module fin_grn_ext_subs

    use, intrinsic :: iso_fortran_env, only: rk => real64
    implicit none

    real(rk), parameter :: pi = acos(-1.0_rk)
    real(rk), parameter :: sqpi = sqrt(pi)

contains

!        purpose: this program computes the error function
!                 erf(x) using subroutine error
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        input:   x   --- argument of erf(x)
!        output:  err --- erf(x)
!        example:
!                   x         erf(x)
!                 ---------------------
!                  1.0       .84270079
!                  2.0       .99532227
!                  3.0       .99997791
!                  4.0       .99999998
!                  5.0      1.00000000

    pure subroutine error(x, err)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: err
        integer k
        real(rk) x2, er, r, c0
        real(rk), parameter :: eps = 1.0e-15_rk

        x2 = x*x
        if (dabs(x) < 3.5_rk) then
            er = 1.0_rk
            r = 1.0_rk
            do k = 1, 50
                r = r*x2/(k + 0.5_rk)
                er = er + r
                if (dabs(r) <= dabs(er)*eps) exit
            end do
            c0 = 2.0_rk*x*dexp(-x2)/sqpi
            err = c0*er
        else
            er = 1.0_rk
            r = 1.0_rk
            do k = 1, 12
                r = -r*(k - 0.5_rk)/x2
                er = er + r
            end do
            c0 = dexp(-x2)/(dabs(x)*sqpi)
            err = 1.0_rk - c0*er
            if (x < 0.0_rk) err = -err
        end if

    end subroutine

!        purpose: this program computes the exponential integral
!                 e1(x) using subroutine e1xa
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        input :  x  --- argument of e1(x)  ( x > 0 )
!        output:  e1 --- e1(x)
!        example:
!                   x        e1(x)
!                 ----------------------
!                  0.0     .1000000+301
!                  1.0     .2193839e+00
!                  2.0     .4890051e-01
!                  3.0     .1304838e-01
!                  4.0     .3779352e-02
!                  5.0     .1148296e-02
!
    pure subroutine e1xa(x, e1)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: e1
        real(rk) es1, es2

        if (x == 0.0_rk) then
            e1 = 1.0d+300
        else if (x <= 1.0_rk) then
            e1 = -dlog(x) + ((((1.07857e-3_rk*x - 9.76004e-3_rk)*x + 5.519968e-2_rk)*x - 0.24991055_rk)*x &
                             + 0.99999193_rk)*x - 0.57721566_rk
        else
            es1 = (((x + 8.5733287401_rk)*x + 18.059016973_rk)*x + 8.6347608925_rk)*x + 0.2677737343_rk
            es2 = (((x + 9.5733223454_rk)*x + 25.6329561486_rk)*x + 21.0996530827_rk)*x + 3.9584969228_rk
            e1 = dexp(-x)*es1/(es2*x)
        end if

    end subroutine

!        purpose: this program computes the exponential integral
!                 en(x) using subroutine enxb
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        input :  x --- argument of en(x)
!                 n --- order of en(x)  (n = 0,1,2,...)
!        output:  en --- en(x)
!        example: x = 10.0
!                    n         en(x)
!                  ----------------------
!                    0     .45399930d-05
!                    1     .41569689d-05
!                    2     .38302405d-05
!                    3     .35487626d-05
!                    4     .33041014d-05
!                    5     .30897289d-05
!
    pure subroutine enxa(n, x, en)
        integer, intent(in) :: n
        real(rk), intent(in) :: x
        real(rk), intent(out) :: en
        integer k
        real(rk) e0, e1

        e0 = dexp(-x)/x

        call e1xa(x, e1)

        if (n == 0) then
            en = e0
        elseif (n == 1) then
            en = e1
        else
            do k = 2, n
                en = (dexp(-x) - x*e1)/(k - 1.0_rk)
                e1 = en
            end do
        end if

    end subroutine

!        purpose: this program computes the exponential integral
!                 ei(x) using subroutine eix
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        input :  x  --- argument of ei(x)
!        output:  ei --- ei(x) ( x > 0 )
!        example:
!                   x        ei(x)
!                 -----------------------
!                   0    -.10000000+301
!                   1     .18951178e+01
!                   2     .49542344e+01
!                   3     .99338326e+01
!                   4     .19630874e+02
!                   5     .40185275e+02

    pure subroutine eix(x, ei)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: ei
        integer k
        real(rk) r, ga

        if (x == 0.0_rk) then
            ei = -1.0d+300
        else if (x <= 40.0_rk) then
            ei = 1.0_rk
            r = 1.0_rk
            do k = 1, 100
                r = r*k*x/(k + 1.0_rk)**2
                ei = ei + r
                if (dabs(r/ei) <= 1.0d-15) exit
            end do
            ga = 0.5772156649015328_rk
            ei = ga + dlog(x) + x*ei
        else
            ei = 1.0_rk
            r = 1.0_rk
            do k = 1, 20
                r = r*k/x
                ei = ei + r
            end do
            ei = dexp(x)*ei/x
        end if

    end subroutine

!        purpose: this program computes the gamma function
!                 gamma(x) for x > 0 or its logarithm using
!                 subroutine lgama
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        purpose: compute gamma function gamma(x) or ln[gamma(x)]
!        input:   x  --- argument of gamma(x) ( x > 0 )
!                 kf --- function code
!                        kf=1 for gamma(x); kf=0 for ln[gamma(x)]
!        output:  gl --- gamma(x) or ln[gamma(x)]
!        examples:
!                   x           gamma(x)
!                 -------------------------
!                  0.5     .1772453851d+01
!                  2.5     .1329340388d+01
!                  5.0     .2400000000d+02
!                  7.5     .1871254306d+04
!                 10.0     .3628800000d+06

    pure subroutine lgama(kf, x, gl)
        integer, intent(in) :: kf
        real(rk), intent(in) :: x
        real(rk), intent(out) :: gl
        integer k, n
        real(rk) x0, x2, xp, gl0

        real(rk), parameter :: a(10) = [8.333333333333333d-02, -2.777777777777778d-03, &
                                        7.936507936507937d-04, -5.952380952380952d-04, &
                                        8.417508417508418d-04, -1.917526917526918d-03, &
                                        6.410256410256410d-03, -2.955065359477124d-02, &
                                        1.796443723688307d-01, -1.39243221690590d+00]

        x0 = x
        if (x == 1.0_rk .or. x == 2.0_rk) then
            gl = 0.0_rk
            if (kf == 1) gl = dexp(gl)
        else if (x <= 7.0_rk) then
            n = int(7 - x)
            x0 = x + n
            x2 = 1.0_rk/(x0*x0)
            xp = 6.283185307179586477_rk
            gl0 = a(10)
            do k = 9, 1, -1
                gl0 = gl0*x2 + a(k)
            end do
            gl = gl0/x0 + 0.5_rk*dlog(xp) + (x0 - 0.5_rk)*dlog(x0) - x0
            if (x <= 7.0_rk) then
                do k = 1, n
                    gl = gl - dlog(x0 - 1.0_rk)
                    x0 = x0 - 1.0_rk
                end do
            end if
            if (kf == 1) gl = dexp(gl)
        end if

    end subroutine

    !> finds the limit of a series
    pure subroutine limes(n, s, rlimes, i_pade, k_pade, err)

! description:
!   finds the limit of a series in the case where only
!   the first n+1 terms are known.
!
! method:
!   the subroutine operates by applying the epsilon-algorithm
!   to the sequence of partial sums of a seris supplied on input.
!   for desciption of the algorithm, please see:
!
!   [1] t. mishonov and e. penev, int. j. mod. phys. b 14, 3831 (2000)
!
! owners: todor mishonov & evgeni penev
!
! history:
! version   date         comment
! =======   ====         =======
! 1.0       01/04/2000   original code. t. mishonov & e. penev
!
! code description:
!   language:           fortran 90.
!   software standards: "european standards for writing and
!                        documenting exchangeable fortran 90 code".
!

!  scalar arguments with intent(in):
        integer, intent(in) :: n            ! width of the epsilon-table

!  array arguments with intent(inout):
        real(rk), intent(inout) :: s(0:n)   ! sequential row of the epsilon-table

!  scalar arguments with intent(out)
        real(rk), intent(out) :: rlimes     ! value of the series limes
        integer, intent(out) :: i_pade      ! power of the numerator
        integer, intent(out) :: k_pade      ! power of the denominator
        real(rk), intent(out) :: err        ! empirical error

!  local parameters                         ! these two need no description ;-)
        real(rk), parameter :: zero = 0.0_rk
        real(rk), parameter :: one = 1.0_rk

!  local scalars
        real(rk) :: a_max       ! maximum element of a
        integer :: i            ! index variable for columns
        integer :: k            ! index variable for rows

!  local arrays
        real(rk) :: a(0:n)      ! auxiliary row of the epsilon-table

!- end of header --------------------------------------------------------------

!  parse input: the algorithm cannot employ more elements than supplied on
!               input, i.e. n <= size(s)
!
        ! if (n > size(s(:))) then
        !     write (*, '(a)') '*** illegal input to limes: n > size(s)'
        !     stop 1
        ! end if

!  algorithm not applicable for n < 2
!
        ! if (n < 2) then
        !     write (*, '(a)') '*** illegal input to limes: n < 2'
        !     stop 2
        ! end if

!  i. initialize with natural assignments

        rlimes = s(n)                       ! the n-th partial sum
        err = abs(s(n) - s(n - 1))          ! error -> |s(n) - s(n-1)|
        i_pade = n                          ! pade approximant [n/0]
        k_pade = 0                          !
        a(:) = zero                         ! auxiliary row initially set to zero
        a_max = zero                        ! max. element set to zero
        k = 1                               ! algorithm starts from the first row

! ii. main loop: fill in the epsilon table, check for convergence ...
!     (provision against devision by zero employs pseudo-inverse numbers)
        do
            if (n - 2*k + 1 < 0) exit

! update the auxiliary row a(i) of the epsilon-table
! by applying the "cross rule".
!
            do i = 0, n - 2*k + 1
                if (s(i + 1) /= s(i)) then
                    a(i) = a(i + 1) + one/(s(i + 1) - s(i))
                else
                    a(i) = a(i + 1)
                end if
            end do
            if (n - 2*k < 0) exit

!  update the sequential row s(i) of the epsilon-table
!  by applying the "cross rule".
!
            do i = 0, n - 2*k
                if (a(i + 1) /= a(i)) then
                    s(i) = s(i + 1) + one/(a(i + 1) - a(i))
                else
                    s(i) = s(i + 1)
                end if

!  check for convergence, based on a_max; see ref. [1]
!
                if (abs(a(i)) > a_max) then
                    a_max = abs(a(i))
                    rlimes = s(i)
                    k_pade = k
                    i_pade = i + k_pade
                    err = one/a_max
                    if (s(i + 1) == s(i)) return
                end if
            end do
            k = k + 1      ! increment row index
        end do

    end subroutine limes

!        purpose: compute modified bessel functions i0(x), i1(1),
!                 k0(x) and k1(x), and their derivatives
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        input :  x   --- argument ( x � 0 )
!        output:  bk0 --- k0(x)
!                 bk1 --- k1(x)
!

    pure subroutine ik01a(x, bk0, bk1)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: bk0, bk1
        integer k, k0
        real(rk) el, x2, bi0, bi1, r, xr, ca, ct, w0, ww, cb, xr2
        real(rk), parameter :: a(12) = [0.125_rk, 7.03125e-2_rk, &
                                        7.32421875e-2_rk, 1.1215209960938d-1, &
                                        2.2710800170898d-1, 5.7250142097473d-1, &
                                        1.7277275025845_rk, 6.0740420012735_rk, &
                                        2.4380529699556e01_rk, 1.1001714026925e02_rk, &
                                        5.5133589612202e02_rk, 3.0380905109224e03_rk]
        real(rk), parameter :: b(12) = [-0.375_rk, -1.171875d-1, &
                                        -1.025390625d-1, -1.4419555664063d-1, &
                                        -2.7757644653320d-1, -6.7659258842468d-1, &
                                        -1.9935317337513_rk, -6.8839142681099_rk, &
                                        -2.7248827311269e01_rk, -1.2159789187654e02_rk, &
                                        -6.0384407670507e02_rk, -3.3022722944809e03_rk]
        real(rk), parameter :: a1(8) = [0.125_rk, 0.2109375_rk, &
                                        1.0986328125_rk, 1.1775970458984e01_rk, &
                                        2.1461706161499e02_rk, 5.9511522710323e03_rk, &
                                        2.3347645606175e05_rk, 1.2312234987631e07_rk]

        el = 0.5772156649015329_rk
        x2 = x*x
        if (x == 0.0_rk) then
            bi0 = 1.0_rk
            bi1 = 0.0_rk
            bk0 = 1.0d+300
            bk1 = 1.0d+300
            return
        else if (x <= 18.0_rk) then
            bi0 = 1.0_rk
            r = 1.0_rk
            do k = 1, 50
                r = 0.25_rk*r*x2/(k*k)
                bi0 = bi0 + r
                if (dabs(r/bi0) < 1.0d-15) exit
            end do
            bi1 = 1.0_rk
            r = 1.0_rk
            do k = 1, 50
                r = 0.25_rk*r*x2/(k*(k + 1))
                bi1 = bi1 + r
                if (dabs(r/bi1) < 1.0d-15) exit
            end do
            bi1 = 0.5_rk*x*bi1
        else
            k0 = 12
            if (x .ge. 35.0) k0 = 9
            if (x .ge. 50.0) k0 = 7
            ca = dexp(x)/dsqrt(2.0_rk*pi*x)
            bi0 = 1.0_rk
            xr = 1.0_rk/x
            do k = 1, k0
                bi0 = bi0 + a(k)*xr**k
            end do
            bi0 = ca*bi0
            bi1 = 1.0_rk
            do k = 1, k0
                bi1 = bi1 + b(k)*xr**k
            end do
            bi1 = ca*bi1
        end if
        if (x <= 9.0_rk) then
            ct = -(dlog(x/2.0_rk) + el)
            bk0 = 0.0_rk
            w0 = 0.0_rk
            r = 1.0_rk
            do k = 1, 50
                w0 = w0 + 1.0_rk/k
                r = 0.25_rk*r*x2/(k*k)
                bk0 = bk0 + r*(w0 + ct)
                if (dabs((bk0 - ww)/bk0) < 1.0d-15) exit
                ww = bk0
            end do
            bk0 = bk0 + ct
        else
            cb = 0.5_rk/x
            xr2 = 1.0_rk/x2
            bk0 = 1.0_rk
            do k = 1, 8
                bk0 = bk0 + a1(k)*xr2**k
            end do
            bk0 = cb*bk0/bi0
        end if
        bk1 = (1.0_rk/x - bi1*bk0)/bi0

    end subroutine

!        purpose: compute bessel functions j0(x), j1(x), y0(x),
!                 y1(x), and their derivatives
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   parameters:
!        input :  x   --- argument of jn(x) & yn(x) ( x � 0 )
!        output:  bj0 --- j0(x)
!                 bj1 --- j1(x)
!                 by0 --- y0(x)
!                 by1 --- y1(x)
!
    pure subroutine jy01b(x, bj0, bj1, by0, by1)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: bj0, bj1, by0, by1
        real(rk) t, t2, a0, p0, q0, ta0, p1, q1, ta1

        if (x == 0.0_rk) then
            bj0 = 1.0_rk
            bj1 = 0.0_rk
            by0 = -1.0d+300
            by1 = -1.0d+300
            return
        else if (x <= 4.0_rk) then
            t = x/4.0_rk
            t2 = t*t
            bj0 = ((((((-.5014415e-3_rk*t2 + 0.76771853e-2_rk)*t2 - 0.0709253492_rk)*t2 + 0.4443584263_rk) &
                     *t2 - 1.7777560599_rk)*t2 + 3.9999973021_rk)*t2 - 3.9999998721_rk)*t2 + 1.0_rk
            bj1 = t*(((((((-.1289769e-3_rk*t2 + 0.22069155e-2_rk)*t2 - 0.0236616773_rk)*t2 + 0.1777582922_rk) &
                        *t2 - 0.8888839649_rk)*t2 + 2.6666660544_rk)*t2 - 3.9999999710_rk)*t2 + 1.9999999998_rk)
            by0 = (((((((-.567433e-4_rk*t2 + 0.859977e-3_rk)*t2 - 0.94855882e-2_rk)*t2 + 0.0772975809_rk)*t2 &
                      - 0.4261737419_rk)*t2 + 1.4216421221_rk)*t2 - 2.3498519931_rk)*t2 &
                   + 1.0766115157)*t2 + 0.3674669052_rk
            by0 = 2.0_rk*dlog(x/2.0_rk)*bj0/pi + by0
            by1 = ((((((((.6535773e-3_rk*t2 - 0.0108175626_rk)*t2 + 0.107657606_rk)*t2 - 0.7268945577_rk)*t2 &
                       + 3.1261399273_rk)*t2 - 7.3980241381_rk)*t2 + 6.8529236342_rk)*t2 &
                    + 0.3932562018_rk)*t2 - 0.6366197726_rk)/x
            by1 = 2.0_rk*dlog(x/2.0_rk)*bj1/pi + by1
        else
            t = 4.0_rk/x
            t2 = t*t
            a0 = dsqrt(2.0_rk/(pi*x))
            p0 = ((((-.9285d-5*t2 + 0.43506e-4_rk)*t2 - 0.122226e-3_rk)*t2 + 0.434725e-3_rk)*t2 &
                  - 0.4394275e-2_rk)*t2 + 0.999999997_rk
            q0 = t*(((((.8099d-5*t2 - 0.35614e-4_rk)*t2 + 0.85844e-4_rk)*t2 - 0.218024e-3_rk)*t2 &
                     + 0.1144106e-2_rk)*t2 - 0.031249995_rk)
            ta0 = x - 0.25_rk*pi
            bj0 = a0*(p0*dcos(ta0) - q0*dsin(ta0))
            by0 = a0*(p0*dsin(ta0) + q0*dcos(ta0))
            p1 = ((((.10632e-4_rk*t2 - 0.50363e-4_rk)*t2 + 0.145575e-3_rk)*t2 - 0.559487e-3_rk)*t2 &
                  + 0.7323931e-2_rk)*t2 + 1.000000004_rk
            q1 = t*(((((-.9173d-5*t2 + 0.40658e-4_rk)*t2 - 0.99941e-4_rk)*t2 + 0.266891e-3_rk)*t2 &
                     - 0.1601836e-2_rk)*t2 + 0.093749994_rk)
            ta1 = x - 0.75_rk*pi
            bj1 = a0*(p1*dcos(ta1) - q1*dsin(ta1))
            by1 = a0*(p1*dsin(ta1) + q1*dcos(ta1))
        end if

    end subroutine

!        purpose: compute bessel functions j0(x), j1(x), y0(x),
!                 y1(x), and their derivatives
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   description:
!     this is a derived version from the subroutine jy01b.
!
!   parameters:
!        input :  x   --- argument of jn(x) & yn(x) ( x � 0 )
!        output:  bj0 --- j0(x)
!                 bj1 --- j1(x)
!
    pure subroutine jy01bj(x, bj0, bj1)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: bj0, bj1
        real(rk) t, t2, a0, p0, q0, ta0, p1, q1, ta1

        if (x == 0.0_rk) then
            bj0 = 1.0_rk
            bj1 = 0.0_rk
            return
        else if (x <= 4.0_rk) then
            t = x/4.0_rk
            t2 = t*t
            bj0 = ((((((-.5014415e-3_rk*t2 + 0.76771853e-2_rk)*t2 - 0.0709253492_rk)*t2 + 0.4443584263_rk) &
                     *t2 - 1.7777560599_rk)*t2 + 3.9999973021_rk)*t2 - 3.9999998721_rk)*t2 + 1.0_rk
            bj1 = t*(((((((-.1289769e-3_rk*t2 + 0.22069155e-2_rk)*t2 - 0.0236616773_rk)*t2 + 0.1777582922_rk) &
                        *t2 - 0.8888839649_rk)*t2 + 2.6666660544_rk)*t2 - 3.9999999710_rk)*t2 + 1.9999999998_rk)
        else
            t = 4.0_rk/x
            t2 = t*t
            a0 = dsqrt(2.0_rk/(pi*x))
            p0 = ((((-.9285d-5*t2 + 0.43506e-4_rk)*t2 - 0.122226e-3_rk)*t2 + 0.434725e-3_rk)*t2 &
                  - 0.4394275e-2_rk)*t2 + 0.999999997_rk
            q0 = t*(((((.8099d-5*t2 - 0.35614e-4_rk)*t2 + 0.85844e-4_rk)*t2 - 0.218024e-3_rk)*t2 &
                     + 0.1144106e-2_rk)*t2 - 0.031249995_rk)
            ta0 = x - 0.25_rk*pi
            bj0 = a0*(p0*dcos(ta0) - q0*dsin(ta0))
            p1 = ((((.10632e-4_rk*t2 - 0.50363e-4_rk)*t2 + 0.145575e-3_rk)*t2 - 0.559487e-3_rk)*t2 &
                  + 0.7323931e-2_rk)*t2 + 1.000000004_rk
            q1 = t*(((((-.9173d-5*t2 + 0.40658e-4_rk)*t2 - 0.99941e-4_rk)*t2 + 0.266891e-3_rk)*t2 &
                     - 0.1601836e-2_rk)*t2 + 0.093749994_rk)
            ta1 = x - 0.75_rk*pi
            bj1 = a0*(p1*dcos(ta1) - q1*dsin(ta1))
        end if

    end subroutine

!        purpose: compute bessel functions j0(x), j1(x), y0(x),
!                 y1(x), and their derivatives
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   description:
!     this is a derived version from the subroutine jy01b.
!
!   parameters:
!        input :  x   --- argument of jn(x) & yn(x) ( x � 0 )
!        output:  bj0 --- j0(x)
!                 bj1 --- j1(x)
!                 by0 --- y0(x)
!                 by1 --- y1(x)
!
    pure subroutine jy01by0(x, by0)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: by0
        real(rk) bj0, t, t2, a0, p0, q0, ta0

        if (x == 0.0_rk) then
            bj0 = 1.0_rk
            by0 = -1.0d+300
            return
        else if (x <= 4.0_rk) then
            t = x/4.0_rk
            t2 = t*t
            bj0 = ((((((-.5014415e-3_rk*t2 + 0.76771853e-2_rk)*t2 - 0.0709253492_rk)*t2 + 0.4443584263_rk) &
                     *t2 - 1.7777560599_rk)*t2 + 3.9999973021_rk)*t2 - 3.9999998721_rk)*t2 + 1.0_rk
            by0 = (((((((-.567433e-4_rk*t2 + 0.859977e-3_rk)*t2 - 0.94855882e-2_rk)*t2 + 0.0772975809_rk)*t2 &
                      - 0.4261737419_rk)*t2 + 1.4216421221_rk)*t2 - 2.3498519931_rk)*t2 + 1.0766115157)*t2 &
                  + 0.3674669052_rk
            by0 = 2.0_rk*dlog(x/2.0_rk)*bj0/pi + by0
        else
            t = 4.0_rk/x
            t2 = t*t
            a0 = dsqrt(2.0_rk/(pi*x))
            p0 = ((((-.9285d-5*t2 + 0.43506e-4_rk)*t2 - 0.122226e-3_rk)*t2 &
                   + 0.434725e-3_rk)*t2 - 0.4394275e-2_rk)*t2 + 0.999999997_rk
            q0 = t*(((((.8099d-5*t2 - 0.35614e-4_rk)*t2 + 0.85844e-4_rk)*t2 &
                      - 0.218024e-3_rk)*t2 + 0.1144106e-2_rk)*t2 - 0.031249995_rk)
            ta0 = x - 0.25_rk*pi
            by0 = a0*(p0*dsin(ta0) + q0*dcos(ta0))
        end if

    end subroutine

!        purpose: compute bessel functions j0(x), j1(x), y0(x),
!                 y1(x), and their derivatives
!   license:
!
!     this routine is copyrighted by shanjie zhang and jianming jin.  however,
!     they give permission to incorporate this routine into a user program
!     provided that the copyright is acknowledged.
!
!   author:
!
!     shanjie zhang, jianming jin
!
!   reference:
!
!     shanjie zhang, jianming jin,
!     computation of special functions,
!     wiley, 1996,
!     isbn: 0-471-11963-6,
!     lc: qa351.c45.
!
!   modified on:
!     25 july 2017, by yingyi liu
!
!   description:
!     this is a derived version from the subroutine jy01b.
!
!   parameters:
!        input :  x   --- argument of jn(x) & yn(x) ( x � 0 )
!        output:  by1 --- y1(x)
!
    pure subroutine jy01by1(x, by1)
        real(rk), intent(in) :: x
        real(rk), intent(out) :: by1
        real(rk) bj1, t, t2, a0, p1, q1, ta1

        if (x == 0.0_rk) then
            bj1 = 0.0_rk
            by1 = -1.0d+300
            return
        else if (x <= 4.0_rk) then
            t = x/4.0_rk
            t2 = t*t
            bj1 = t*(((((((-.1289769e-3_rk*t2 + 0.22069155e-2_rk)*t2 - 0.0236616773_rk)*t2 + 0.1777582922_rk) &
                        *t2 - 0.8888839649_rk)*t2 + 2.6666660544_rk)*t2 - 3.9999999710_rk)*t2 + 1.9999999998_rk)
            by1 = ((((((((.6535773e-3_rk*t2 - 0.0108175626_rk)*t2 + 0.107657606_rk)*t2 &
                        - 0.7268945577_rk)*t2 + 3.1261399273_rk)*t2 - 7.3980241381_rk)*t2 &
                     + 6.8529236342_rk)*t2 + 0.3932562018_rk)*t2 - 0.6366197726_rk)/x
            by1 = 2.0_rk*dlog(x/2.0_rk)*bj1/pi + by1
        else
            t = 4.0_rk/x
            t2 = t*t
            a0 = dsqrt(2.0_rk/(pi*x))
            p1 = ((((0.10632e-4_rk*t2 - 0.50363e-4_rk)*t2 + 0.145575e-3_rk)*t2 - 0.559487e-3_rk)*t2 &
                  + 0.7323931e-2_rk)*t2 + 1.000000004_rk
            q1 = t*(((((-0.9173d-5*t2 + 0.40658e-4_rk)*t2 - 0.99941e-4_rk)*t2 + 0.266891e-3_rk)*t2 &
                     - 0.1601836e-2_rk)*t2 + 0.093749994_rk)
            ta1 = x - 0.75_rk*pi
            by1 = a0*(p1*dsin(ta1) + q1*dcos(ta1))
        end if

    end subroutine

end module fin_grn_ext_subs
